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, Although density functional theory (DFT) in principle includes even long-range interactions, 

^^ , standard implementations employ local or semi-local approximations of the interaction energy 

£^ ■ and fail at describing the van der Waals interactions. We show how to modify a recent density 

functional that includes van der Waals interactions in planar systems [Phys. Rev. Lett. 91, 126402 

j^ , (2003)] to also give an approximate interaction description of planar molecules. As a test case we 

use this modified functional to calculate the binding distance and energy for benzene dimers, with 
the perspective of treating also larger, flat molecules, such as the polycyclic aromatic hydrocarbons 

^ : (PAH), 

o 
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T-H ! 
K*- ■ I. INTRODUCTION 

^ , 

Xf^ _ The benzene dimer is the prototype for aromatic interactions and has been studied extensively 

\0 ' both by theoretical Q,Q,y,|j,|^|0| and experimental fT','©,^ means. Describing benzene interactions 

^D ' can be regarded as the first step toward describing interactions of polycyclic aromatic hydrocarbons 

'j . (PAHs). PAHs are planar molecules consisting of several aromatic rings, where the peripheral carbon 

atoms are covalently bonded to hydrogen atoms. Both benzene and PAHs structurally resemble 
graphite and exhibit very similar intra- and intermolecular bond lengths, particularly for (stacks of) 
large PAHs. Like the interactions between the sheets in graphite, the interactions between parallel 
benzene molecules or layers in PAH stacks are dominated by the weak, nonlocal van der Waals 
i-rt ' (vdW) interaction. However, the description of the vdW interaction is absent in the traditional 

P5 ' implementations of density functional theory (DFT), implementations which have otherwise been 

O ■ very successful in describing dense, hard materials on the atomic scale. 

J-^ ' The past few years a number of publications [Tfl llil IH, [l3> llJ> HI; H^ H^ uM have addressed 

the problem of consistently extending the common DFT implementations to also include the vdW 
1 J ■ interaction. For the mutual interactions of planar sheets or surfaces a functional was obtained in 

rS I Refs. [i^jll^- For example, a relatively good estimate of the binding distance between two graphite 

planes (graphene) was obtained with this planar functional |l5| , while the usual generalized gradient 
approximations (GGA) fail to bind the graphene layers or do so at an unphysical binding distance and 
energy. Since benzene and PAHs are similar to graphite, the correction to the regular calculations 
will be significant also in dimers of these molecules. 

In this work, we study dimers of benzene (CgHg) by modifying and approximating the planar 
nonlocal density functional [iSllal in order to treat molecules of finite size. We here concentrate on 
dimers where the molecules are placed directly on top of each other ("AA stacking"). The method 
we use is developed with the aim of being able to study also PAH molecules, and it is described 
below. The results for some of the small PAH molecules naphthalene, anthracene, and pyrene 
(CioHg, C14H10, and CigHio) will be reported in a forthcoming publication |ig |. 
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II. METHOD 

DFT in principle does include even long-range interactions, such as the vdW interactions, but 
those are a part of the exchange-correlation functional Er^,, which is approximated in the implementa- 
tions of DFT. The most commonly used approximations to E-j;,, are the local density approximation 
(LDA) and the generalized gradient approximation (GGA), both of which are local or semi- local 
and do not include the vdW interaction. However, GGA (and to some extend LDA) has proven 
very successful in describing the interactions within dense matter, e.g., within a molecule, and a 
consistent correction scheme that includes the long-range intermolecular vdW interactions must not 
change the short-range intramolecular interaction. 



A. The graphite interactions 

In the vdW correction scheme for planar sheets and surfaces, defined and described in Refs. 
|13Lll5lll6| . the exchange energy from GGA is retained, and the correlation energy is determined as 
a sum of the local correlation (obtained from LDA) and the nonlocal correlation £""' , which will be 
described below. For the exchange part the GGA rcvPBE flavor of Zhang and Yang |2^ is chosen 
because it is fitted to exact exchange for atoms. Thus the (vdW-corrected) total energy is written 
as 

-EvdW-DF — Eqqa — Eqqa.c + -EldA.c + E^ , (1) 

where all energy terms are functionals of the electron charge density n(r), and -Bgga.c and i?LDA,c are 
the correlation from GGA and LDA, respectively. These energy terms are directly available from 
our DFT calculations. Our DFT calculations are performed by the plane- wave pseudo-potential 
based program Dacapo |2y . 

The nonlocal correlation, per area A, for a planar, translationally invariant system such as 
graphite is given by [l^ 

where 0(z) fulfills the differential equation (ekc/)')' = ^^£^0 with boundary conditions (/)(0) — and 
(f>{L) — 0, and (j)o{z) is the vacuum solution. The system is enclosed in a long box of length L 
along the direction perpendicular to the graphite plane(s). Cfe is the dielectric function of graphite, 
Fourier-transformed in the plane of the graphite sheet 



+ 4(z)(P + qlf/S + (fc2 + qlf/A 



efc(^, «m) — 1 + __,2 , ,.2r^\n„2 , ^2 \2 /o , (1,2 , „2 \2 /a ■ i^) 



The plasmon frequency Up is defined by i^p{z) = 47rn(z), the Fermi velocity by vf(z) — (37r^n(z))3 , 
n{z) is the planarly averaged electron density varying in the z-direction perpendicular to the plane, 
and iu is the imaginary frequency. We use atomic units unless otherwise noted. 

The parameter g_L is introduced to compensate for the locality of e^ in the z-direction. The 
value of q±^ is materials dependent, and it is found from a separate set of GGA DFT calculations 
with an applied, static electric field across the graphite layer. 

In Figure^a) we show the interaction curve of graphene-graphene in the A A stacking according 
to C3). These results of graphene in the A A stacking allow us to directly compare to our results of 
benzene, and later PAHs, also in AA stacking. The binding separation 4.0 A of graphene (AA) is 
larger than for graphene in the physically correct AB stacking (3.76 A), calculated within the same 
vdW correction scheme |l3- We used the value q±^ = 0.756 in atomic units for a graphene sheet 



B. Extending the functional to treat planar molecules 

It is a natural next step to modify i5"' for planar, translational invariant sheets to treat planar, 
large, but not infinitely extended molecules. The PAH molecules are pieces of graphene, passivated 



by hydrogen atoms at the broken carbon-carbon bonds. We would therefore like to introduce 
appropriate approximations in order to apply the functional to planar, large, finite molecules. 

The (semi-)local part of the total energy in |^, i^cGA — ^'gga.c + ^lda,c, can be determined 
directly from a usual DFT implementation. Such calculations provide us, besides the energy terms, 
with a three-dimensional electron density n(r). In the formalism for £""' in Q the density profile 
n(z) is assumed translationally invariant along the molecule. For the finite molecules we take 
n{z) = ^~oi / dxdy n{x, y, z) where Amoi is the size of the molecular area that affects the vdW 
interactions; this area must be estimated. 

For carbon-dominated molecules such as benzene and PAH the dominating contribution to the 
vdW correction comes from the carbon atoms. We thus assume that the dielectric function for the 
molecules is well described by the functional form ||2Il with the graphene value q±^ — 0.756 a.u. By 
comparing the dacapo-calculated static [u — 0) electric response of the isolated molecule with that 
found from the model dielectric function Q, with Amoi as a parameter, we can determine Amoi- For 
benzene we find that Amoi = 51.7 A^ reproduces the static response. 

Given the model dielectric function ^ with the parameter Amoi determined as described 
above, the calculation of the E"^^ contribution follows from ^. The physical analogy of using the 
molecularly averaged n{z) in ((SJ is that of molecules moved together such that they each take up the 
lateral area Amoi. £^"' for molecules calculated this way therefore somewhat overcounts the vdW 
interaction by including not only the physically correct interaction of one dimer molecule with its 
partner, but also the interaction with the partner's neighbors. Because the vdW attraction falls off 
with distance this effect decreases as the molecule size increases. The other energy terms of (^ are 
not affected. 

We stress that Amoi is larger than the area of the same amount of carbon atoms in a graphene 
plane. In our procedure of finding Amoi we compare the (static) response of a molecule isolated 
in all three directions. We make sure that our Dacapo-calculated response really does come from 
an isolated molecule by checking the convergence of the response with increasing lateral unit cell 
size. The response of the low-density charge distribution region within the molecular plane is thus 
included in our description, and therefore the vdW interaction of these sheets of close molecules is 
different than the vdW interaction of two sheets of graphite, even though the molecules are mostly 
composed of carbon atoms. 

With the assumptions introduced above, necessary for applying the vdW planar functional to 
finite molecules, the quality of the vdW interaction calculations increases with increasing molecular 
size. We emphasize that in this sense, the benzene molecule is an extreme molecule, and we do not 
expect full quantitative agreement of our results with full-powered quantum-chemistry calculations 
mial or with results in our new (more costly) general-geometry formulation of the vdW-DF Iq. 



III. RESULTS 

We here present the results of applying the modified vdW functional to the benzene dimer in 
the AA stacking. The first step is to determine the structure of the isolated molecule, the static 
response to an applied electric field, the (semi-)local energy contributions to the dimer binding, 
and the electron density of the dimer. This step is carried out in a standard DFT program. The 
second step is to calculate the nonlocal correlation energy contribution, i?"' , with the approximations 
discussed above. 

The standard DFT calculations were performed with the plane-wave pseudopotential code 
Dacapo 21] with periodic boundary conditions. We performed the self-consistent calculations (before 
adding S"') using the revPBE GGA functional. The choice of unit cell size involves a trade-off 
between manageable calculations on the one hand, and large lateral and vertical separations between 
periodically repeated images of the molecules (or dimers) on the other hand. When deciding on the 
unit cell size we also benefit from the short-range nature of the traditional DFT-implementations. 
In our calculations we used a hexagonal unit cell of size (17.112 A, 17.112 A, 26 A). 

In the plane-wave code the number of plane waves is determined by a cutoff in the plane- 
wave energy. A preliminary analysis showed a cutoff at 450 eV to be sufficient. However, as we 
discuss below, there is a small energy difference between dimers within the same unit cell at very 
large separation and two isolated molecules in two separate unit cells. This offset depends on the 
plane-wave cutoff in a non-trivial way. We emphasize that the offset in the total energy difference 
is unrelated to i?"'. 
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FIG. 1: Energy curves for graphene and benzene dimers in the AA stacking as a function of separation dis- 
tance. Panel (a) shows the graphene-graphene interaction, panel (b) shows the benzene-benzene interaction. 
The dashed lines are the self-consistently determined revPBE total energies, the dotted lines show the i?"' 
correction, and the solid lines show the _Evdw-DF total energy curves obtained according to Eq. Q. The 
circles on the E^^w-df curves indicate the calculations performed. 



A. Benzene-benzene interaction energy 

The benzene dimer interaction is calculated by fixing the relative atomic positions within of 
each of the two molecules of the dimer and varying the distance between them, while keeping the 
size of the unit cell constant. Figure ^b) shows the total-energy curves for benzene: the revPBE 
curve (which shows no binding at all), the nonlocal correlation E^\ and the full i^vdW-DF that 
includes nonlocal correlation. The reference energy is the benzene dimer at large (11 A) separation. 

We find in vdW-DF the binding energy 100 meV at the distance 4.1 A. These results are in 
agreement with our results for the AA stacked graphene dimer as we expect the binding distance 
there to be similar to that of the benzene dimer and PAH dimers in the AA stacking. Both the 
benzene separation distance and interaction energy are in reasonable, although not perfect agreement 
with experiments and other calculations. For instance, Refs. Q, |j, |^ report calculating binding 
separations of around 3.8 A and binding energies in the range of 64-92 meV. As mentioned earlier, 
we expect the quality of our results to improve for larger molecules. 

As previously mentioned, no binding arises in the revPBE flavor of GGA. We note that GGA 
Perdew-Wang 91 (PW91) 22] calculations (not shown) exhibit an unphysically small binding energy 
of the benzene dimer (12 meV) at separation 4.7 A. However, this PW91 binding arises from the 
unphysical interactions mediated by the exchange contribution in the manner previously documented 
for graphite |l5l |. 



B. Discussion of the energy reference level 

During the course of this benzene and PAH study we encountered a technical problem in 
determining an appropriate reference level for our (traditional) total energy calculations based on 
GGA. 

All dimer energies presented in this paper, whether revPBE or vdW-DF energies, are given 
with reference to the energy of two molecules far apart. Thus the energies that we report, and which 
we for now generally denote E'diff , are total energy differences given by 

E^isid) = Etotid) - E,,i (4) 

where i^tot (d) is the total energy of the dimer determined at separation d, and Ej-^i is the energy 
of two molecules far apart. This applies both to the benzene molecules and the graphite sheets. 
The reference level E^^i can be calculated either by removing one molecule to find half the reference 
energy, or by keeping both molecules within the unit cell and moving them "far apart" within the 
box. In both cases, we make sure that the periodic cell is sufficiently large that the molecule does 
not interact (electrostatically) with its periodic images. In our underlying GGA-DFT calculations 
it turned out that these two seemingly equivalent ways of determining E^^i lead to two slightly 
different results. Numerically, the difference is more significant for the benzene molecules than for 
graphene. 

In Fig.^b) we report curves of interaction energies of dimers with increasing separation. The 
natural choice of reference energy i^iof is therefore the dimer with the two molecules "far apart" , by 
which we mean 11-12 A separation in boxes of length of 24-26 A. However, the problem merits an 
analysis. 

In Figure [3 we show the difference in reference energies SE,-ci — defined as the energy of a 
pair of molecules "far apart" subtracted by the energy of two separate molecules — as a function 
of the choice of DFT convergence parameters. In Dacapo, besides the plane-wave energy cutoff for 
the wave functions, it is possible to set a different (larger) cutoff for the charge density, to enhance 
the description of the charge density at the price of only a modest increase in the calculational size. 
Most of the calculations in Figure |21 use a larger charge-density cutoff than the wave- function cutoff. 
In order to analyze SE^et also for large charge-density cutoff energies we carried out a number of 
calculations in a unit cell with reduced lateral size (10 A). 

A few immediate conclusions can be drawn on the basis of Figure |5] By comparing the single- 
cutoff and the double-cutoff results at charge-density cutoff 500 eV and 550 eV we find that the 
quality of the wave- function description does not affect SE,-ci. We see this in Figure|21by the collapsed 
data points (indicated by arrows) for the revPBE total energy, but it is true also for the exchange 
and the correlation parts of the energy, as well as for PW91 (not shown). This indicates that the 
(noninteracting) kinetic part of the total energy, which is calculated from the wave functions, does 
not contribute to dE^-ef. 

In contrast, SEy-ef depends non-trivially on the charge-density cutoff. It is not clear that a very 
high cutoff value will yield good results, rather, SE^ci seems to oscillate with charge-density cutoff. 
By comparing the results for the usual box size with the small-box results we further notice that 
(5-Ercf decreases with the amount of low-density (vacuum) region in the box. SE^ci thus seems to be 
mainly affected by contributions from the low-density regions. A major part of SEj-d comes from 
the exchange part of the total energy. Not shown in the figure are SEj-d for the PW91 total energy 
and for correlation (revPBE and PW91). These energies all give a finite, but smaller value of SE^ci- 
Unfortunately the problems are worse exactly for the revPBE-GGA that enters our subsequent 
vdW-DF calculations in (^j. 

Besides the ever-present risk of undetected programming errors, we speculate on three pos- 
sible origins for the finite SE^^i that affects our underlying GGA-DFT calculations: (i) numerical 
instabilities and formal inaccuracies arising from the implementation of the exchange-correlation 
approximation; (ii) the mixing of different exchange-correlation approximations for creation of pseu- 
dopotentials and use in calculations; and (iii) the use and/or implementation of ultrasoft pseudopo- 
tentials. The possibility of formal problems in the exchange-correlation term borrows some support 
from the work by Lacks and Gordon in 1993 z3]: PW91 exchange, and a number of other exchange 
approximations, were shown to give values of |(5_Ei.cf | for noble gas dimers far exceeding those from 
exact exchange. The deviation arises from contributions in the low-density, large-gradient regions 
of space, which is precisely the regions that dominate our system. If this will turn out to be the 
origin ^^. 




400 600 800 1000 1200 1400 1600 1800 2000 
Charge-density cutoff energy [eV] 



FIG. 2: DifTerence in reference energies (5-Erof as a function of charge-density cutoff energy. Diamonds are for 
calculations with wave-function cutoff energy equal to the charge-density cutoff energy, empty circles have 
wave-function cutoff 450 eV in the usual unit cell size and black circles have wave-function cutoff 450 eV in 
a laterally smaller unit cell, as described in the text. Solid, dashed, and dash-dotted lines show 5-Bref of the 
revPBE total energy, of revPBE exchange, and of PW91 exchange. 



of at least a part of the offset SE^ai, it is not too surprising that this does not seem to be 
perceived as a general problem, since GGA and most present-day DFT codes have been implemented 
and tested mostly on dense systems with less vacuum. 

We might worry that using a mix of exchange-correlation approximations in the pseudo- 
potential and in the energy calculations could lead to unwanted effects. Our calculations within 
traditional DFT were carried out with Vanderbilt ultrasoft pseudopotentials [23, originally cre- 
ated within the PW91 flavor of the GGA functional, and the energy and the charge density were 
self-consistently determined within revPBE. However, we carried out a set of calculations for the 
graphene-graphene system, where the offset is much less in size, that showed this mixing to not be 
the cause of the offset. The calculations of SEj-d were carried out with all combinations of pseudopo- 
tentials (revPBE, PW9I, PBE ^), self-consistently determined charge density (revPBE, PW91, 
PBE), and energy extracted (revPBE, PW91, PBE, RPBE ^, LDA, as well as all exchange and 
correlation contributions). In all cases, the major part of the offset arises from the exchange part 
of the energy, whereas the offset on the GGA and LDA correlation is generally almost an order of 
magnitude smaller. 

We cannot, within this study, test whether the use of ultrasoft pseudopotentials contributes 
to 6Ere{- This issue requires controlled comparison of results from different codes and choice of 
pseudopotentials and will be addressed in a forthcoming investigation. 

In addition to the tests mentioned above, we also tested that the exact position of the monomer 
with respect to the Fast-Fourier- Transform grid does not matter on the scale of energy differences 
considered here. Further, we speculated that SE^d might be due to numerical errors arising when 
very small charge-density values are left out in the calculation of the energies, but by changing the 
limit of the smallest charge-density values included we found this to be irrelevant to 5E,-c{- 



IV. CONCLUSIONS 



We have presented a modification of an existing planar-system van-der-Waals DFT functional 
to approximately treat also large but finite planar molecules. We use the modified functional for the 
benzene dimer as an example, with the prospect of applying it to polycyclic aromatic hydrocarbon 



7 

molecules. We have characterized a problem in traditional implementations of DFT for molecular 
systems with large volumes of low charge density. This problem motivates further investigations. 
The modified vdW functional presented is a step towards being able to treat sparse matter systems 
consistently within DFT. 
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